Polymorphism, thermodynamic anomaUes and network formation in an atomistic 

model with two internal states 
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Using molecular dynamics simulations we study the temperature-density phase diagram of a 
simple model system of particles in two dimensions. In addition to translational degrees of freedom, 
each particle has two internal states and interacts with a modified Lennard- Jones potential which 
depends on relative positions as well as the internal states. We find that, despite its simplicity, the 
model has a rich phase diagram showing many features of common network-forming liquids such 
as water and silica, including polymorphism and thermodynamic anomalies. We believe our model 
may be useful for studies concerning generic features of such complex liquids. 
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I. INTRODUCTION 

Water pi, molten silicon [2] and silica glasses [3 are com- 
mon examples of network-forming liquids which are dis- 
tinguished from simple fluids (like liquid argon) by the 
presence of complex intermolecular interactions allow- 
ing bond-formation only in favoured directions. Simi- 
lar ability is also seen in macro-molecular systems such 
as colloidal particles with patchy interactions [H |5]. Due 
to their propensity toward forming extensive networked 
structures, especially in supercooled states, such liquids 
show many intriguing properties such as isobaric density 
maxima, negative thermal expansion coefficients (ofp), 
and increase in the isothermal compressibility (xt) and 
isobaric heat capacity(cp) upon supercooling. Also, 
network- formers exhibit polymorphism i.e. existence of 
multiple crystalline as well as amorphous states varying 
in local coordination and density [6 . Over the years a 
lot of theoretical, computational and experimental ef- 
fort has been made in order to understand the origin 
of such anomalous properties of network- formers espe- 
cially water |7]-[T2]. It is now well established from a mi- 
croscopic point of view that the anomalous properties of 
water are related to the existence of high and low den- 
sity forms [1, 13, JA^ of the supercooled liquid in addi- 
tion to the several stable states of ice. The met ast able 
liquid-liquid critical point [15, 16 in water, buried deep 
within the (stable) ice phase can, nevertheless, crucially 
influence the dynamics of nucleation and amorphization. 
A very similar scenario has recently been proposed for 
liquid silicon [17 which suggests that thermodynamical 
anomalies may be expected for any substance which has 
low, as well as high, density crystalline states. Since low 
energy open crystals are most easily (though not exclu- 
sively) formed in systems of molecules with highly direc- 
tional bonding e.g. the hydrogen bond network in water, 
strongly orientation dependent interactions appear to be 
a prerequisite for water-like properties [18 . 

While the connection between directional bonding, 
polymorphism, liquid-liquid phase transitions and ther- 
modynamic anomalies is established, the univeral na- 
ture of this connection opens up many new avenues 
of investigation. What is the effect of confining fields 
and substrates on the properties of networked liquids? 



What happens to the network driven under shear flow? 
How does one characterize the dynamics of network for- 
mation? To answer these, and other similar questions 
one ideally needs a simple and generic model which is 
amenable to simulation on the computer, as well as acces- 
sible to fairly simple theoretical analyses. Though such 
generic models do exist, they are either idealized, as in a 
lattice model[18, 19 , or not much simpler to handle than 
fully realistic water (or silica) models [20]; the most time 
consuming part being an accurate treatment of molec- 
ular rotations and translation-rotation coupling at high 
densities [2T]. 

In this paper, we describe a simple and generic atom- 
istic model system which shows many of the characteris- 
tics of network forming liquids without the complication 
of molecular rotations. Our model consists of particles 
i = 1 . . . N which have an internal coordinate Si . Unlike 
an angular coordinate. Si are discrete and can take only 
one of two possible values ±1. The interaction between 
particle pairs z, j is strongly directional and depends on 
the displacement Vij as well as the values of Si and Sj. 
The internal state Si may be thought of as mimicking 
rotations though it is not necessary to do so. 

Our main results are the following. Despite the 
simplicity of the model, we obtain a rich phase di- 
agram (Figjll) using finite-size scaled block analysis 
technique [2 2j. The system has two distinct crystalline 
phases viz. a low density honeycomb lattice (HS) and 
a closed packed triangular crystal (TS). These crystals 
melt into a liquid (L) or sublimate to a vapor (G) under 
appropriate conditions (FigJ2[a)-(d)). Remarkably, the 
HS phase coexists with a liquid which is of higher den- 
sity than the solid; very similar to water and (normal) ice 
at atmospheric temperatures and pressures. The liquid 
phase at low temperatures shows anomalous variation of 
pressure and density with temperature similar to that 
of water and silicon. Correlations in the liquid phase |23] 
begin to show strong directional modulations before such 
anomalies arise pointing to the formation of a dynamic 
network within the liquid. Because of the simplicity of 
our model, such correlations can easily be obtained from 
perturbative, liquid state integral equation theory. 

This paper is organized as follows. In Section 2 we 
introduce the essential features of the model and provide 



details of our simulations. We discuss our results focusing 
on the complete density-temperature phase-diagram of 
the system, anomalous properties and the development 
of short-range order in Section 3. We conclude in Section 
4 pointing out some directions of future work. 
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FIG. 1: (color-online) Complete phase diagram of the system 
in the p — T plane as obtained from our simulations. The 
red circles, obtained from block- analysis technique, show the 
gas (G) and hexagonal-solid (HS) densities at different tem- 
peratures. The black squares, obtained from pressure-density 
isotherms, show the phase boundaries of HS-liquid (L), HS- 
triangular-solid(TS) and L-TS coexistance regions. The solid 
lines joining the red circles and black squares are guides to the 
eye. The two phase coexistence regions are shaded in gray. 
Also shown is the locus of the pressure-minima as a dotted 
line. 




FIG. 2: Snapshot picture after 10^ MD-steps for (a) p — 
0.3, T = 0.1 showing the G-HS co-existence, (b) p — 0.6, T — 
0.1 showing the HS phase, (c) p = 0.7, T = 0.3 showing the L 
phase, (d) p — 0.86, T = 0.1 showing HS-TS co-existance. 



II. THE MODEL AND SIMULATION DETAILS 

A. Interaction potential 

Our system consists of N particles in a box of area A, 
with 50 : 50 mixture of particles with the internal co- 
ordinates Si = ±.1 interacting with each other through 
a simple pairwise-additive 2-body interaction in two di- 
mensions. The interaction potential between particles i 
and j, separated by the radius vector r^j, viz. UsiSji'^ij) 
may be decomposed as follows: 
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where r = |r^j| is the magnitude and = cos~^ (vij-i/rij) 
the angle r^j makes with the x-axis. The isotropic {0 
independent) parts of the potential are given by, 



Voir) = ULj{r)^e, r<r^ 



(2) 



with ULj{r) = 4e[(cr/r)i2 - (a/r)^], the usual 6-12 
Lennard- Jones interaction; r^in, = 2'^a and e is the po- 
sition and depth of the minimum in Ulj respectively and, 



where, 



Vi{r) = Vo{r)^Vi{r) 



2Vi{r) = -e-Ecut, r < rr, 



ULj{r) - Ecut, r > rr, 
angle dependent part, 



(3) 



(4) 



with Ecut = ULj{rcut) with r^„j = 2.5(7^. Finally, the 



SV{r,e) = Vi{r)cos{W) 



(5) 



In FigJ3]we have plotted the full interaction potential 
as a function of r^j for two cases, namely when Si = Sj = 
1 (red/light grey curve) and Si = —Sj = 1 (blue/dark 
grey curve) along the direction 6ij = 0. The form of 
the potential ensures that an interchange of the labels i 
and j is accompanied by a rotation of the potential by tt. 
The units for distances and energies are the same as in 
a Lennard- Jones system, viz. a and e and for simplicity 
we have taken <j = e = 1 without loss of generality. 



B. Details of the simulations 

We perform molecular dynamics simulation [21] with 
N = 1176 particles in a two-dimensional rectangular box 
of area A using a standard velocity -verlet algorithm tak- 
ing care to choose an integration time step At = 10~^ 
such that the total energy(E) is conserved to 1 in 10^ 




FIG. 3: (color-online) Plot of the interaction potential be- 
tween two particles with same: red/light grey curve (1,1), 
and opposite: blue/dark grey curve (1,-1) values for {Si, Sj). 
The plots are along the -\-ve x-axis where the potential has 
the deepest minimum. Inset top left: Polar plot of the an- 
gular part of the potential for Si = 1 and Sj = —1 showing 
the three directions along which the potential has minima. 
Inset top right: the geometries for the two distinct cases are 
pictorially illustrated. 



or better [2T. We perform equilibration runs in constant 
N,A,T ensemble, while production runs are performed in 
the constant N,A,E ensemble. The temperature is fixed by 
rescaling the velocities and thermodynamic quantities are 
obtained after equilibration starting from either a hon- 
eycomb or a triangular crystal depending on the density. 
Periodic boundary conditions are applied in both direc- 
tions. In each state, we equilibrated the system for 6x 10^ 
time steps and then calculated time-averages for 4 x 10^ 
time steps. 

We obtain most of the phase diagram shown in Figjl] 
using the block analysis technique [22 which has been 
used in the past for similar systems and has been de- 
scribed in detail. Briefly, the simulation box is divided 
into M5 equal sized blocks and the local density is evalu- 
ated for each block. From the block densities one can con- 
struct the block density probability distribution Pubip) 
for every block indexed by M5. In a single phase region, 
this distribution function is a simple Gaussian centered 
about the mean density p with a width which depends on 
M5 and the isothermal compressibility of the system xt- 
In a co-existence region, the distribution function consists 
of a sum of two Gaussians each centered on the densi- 
ties of the co-existing phases, which are then obtained 
by fitting a sum of appropriately chosen Gaussians and 
extrapolating the resulting co-existence densities to the 
thermodynamic limit M5 ^0. In order that the block- 
analysis technique works, the co-existing densities need 
to be well separated. The block analysis technique there- 
fore works best for liquid -gas and solid -gas coexistence. 
For parts of the phase diagram where the co-existence 
regions are narrow, we obtain co-existing densities in 



the usual way from pressure isotherms using a Maxwell's 
equal-area construction. Apart from co-existence densi- 
ties, we also calculate pair distribution functions for like 
and unlike Si to characterize the various phases. 
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FIG. 4: (color-online) The probability distribution curves for 
density PM^ip) for a few values of Mb in the (a) in the HS 
region(p = 0.628, T = 0.12) and (b) G-HS co-existence region 
{p — 0.35, T = 0.15). The lines joining the points are guides 
to the eye. In (b), we have included two Gaussians (dashed 
red curves) which were fitted to the data for Mb = 14 



III. RESULTS AND DISCUSSIONS 

A. Phase- Diagram 

We obtain the density-temperature phase diagram for 
this system using the block-analysis [22] technique along 
with finite-size scaling described above. In FigH (a) we 
show the density distribution curves for various values of 
M5 in a region where a single homogeneous phase (the 
honeycomb solid, HS) is stable. Each of the curves may 
be represented by single Gaussians with the curves be- 
coming sharper as M5 decreases as expected. Since we 
work in the constant density ensemble, the correspond- 
ing distribution function for M5 = 1 is a trivial delta 
function at p. In contrast, when the system is in a two - 
phase region, we expect the density distribution to be 



bi-modal whenever the size of the blocks are compara- 
ble to the size of the heterogeneities. This is illustrated 
in Figj4] (b) where we have obtained Pubip) fo^" G-HS 
coexistence. For the largest block size shown i.e. for 
M5 = 12, the distribution is still uni- modal, although it 
develops a prominent shoulder. For M5 ^14, however, 
the two peaks can be clearly resolved which gives us the 
coexisting gas and solid densities. 

It is possible that the difference in the coexisting den- 
sities is so small that the bi-modal structure of the den- 
sity distribution cannot be resolved. This happens when 
the coexisting phases are liquid and solid or two differ- 
ent solids. In this case, as mentioned, we use pressure 
isotherms to determine the coexistence densities. The 
pressure isotherms across the L-TS coexistence region 
are shown in Figjs] for three temperatures. The com- 
puted coexistence densities lie along the dashed curves 
as shown. 

The complete density-temperature phase-diagram is 
shown in Figjl] and snapshot pictures of a few of the 
featured phases are shown in Fig|2j Despite the relative 
simplicity of the interactions, our system shows a rather 
complex phase behavior. Firstly, there are two very dif- 
ferent crystalline structures, the honeycomb and trian- 
gular solids, HS and TS. While the former has an open 
structure with strong directional bonding, the latter is 
close packed. The HS solid can coexist with either a low 
density gas, G, or with a liquid, L, which has a density 
higher than itself. This unique feature is reminiscent of 
many networked liquids such as water. The triangular 
solid TS, can be in coexistence with L or even with an- 
other solid viz. HS. Again this is a feature common in 
the water-ice system where various low and high density 
forms of ice may coexist with each other for appropriate 
values of the pressure. 
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B. Thermodynamic anomalies 

Networked liquids which coexist with solids of lower 
density show anomalous behaviour in many thermody- 
namic quantities. To investigate this we plot pressure- 
temperature curves along different iso- chores in FigJoFa) 
and each of these curves show a minimum. A minimum 
in the pressure iso-chore is equivalent to a maximum 
in density along an isobar as can be easily deduced as 
follows [25 . Noting that 



fdP\ _ ap 



(6) 



where xt is the isothermal compressibility, ap is the 
volume-expansion coefficient, it is easy to see that an 
extremum in the isochore is related to vanishing ap. To 
determine the nature of the extremum we consider the 
second derivative. 
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(7) 



Thus, a density maximum (positive dap/dT) leads to 
a pressure minimum. Further, xt itself may show an 
extremum at the point of minimum pressure, a fact which 
is borne out in our system and shown in FigJGJ (b). 
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FIG. 6: (color-online) (a) P{T) isochores at p = 0.66 (red), 
0.68(blue), 0.70(pink) and 0.72(green). Dashed lines are 
drawn as a guides to the eye. (b) The isothermal compressibil- 
ity XT at p = 0.67. Note that the minimum in the pressure as 
well as in the compressibility is very similar to many common 
network forming liquids. 



C. Density correlation functions 



FIG. 5: (color-online) P{p) isotherms along T = 0.2 (red), .25 
(green) and .3(blue). The flat region is the L-TS coexistence 
region. The estimated coexistence densities lie on the dashed 
curves as shown. 



Thermodynamic anomalies arise in a liquid as a re- 
sult of a metastable liquid-liquid phase separation be- 
tween a low density liquid with strong directional corre- 
lations and a more isotropic high density variant at tem- 
peratures where the solid phase is stable ^IG, J/T^. One 



therefore expects that as the temperature is reduced, the 
homogeneous hquid would develop short ranged correla- 
tions which are strongly orientation dependent. In or- 
der to investigate this in our system, we have computed 
density-density correlation functions both from our sim- 
ulations and from a liquid-state integral equation theory. 
In two dimensions, the radial-distribution function ^(r) 




FIG. 7: (color-online) Radial distribution functions gi,i{T) at 
T = 0.5 obtained from (a) simulations and (b) theory and 
5'i,-i(r) at the same temperature obtained from (c) simula- 
tions and (d) theory. 

is defined as the probability of finding a pair of parti- 
cles within r and r + dr of each other given that one of 
particles is at the origin. In our case, we need distribu- 
tion functions for both like, ^1,1 (r), and unlike, ^i^_i(r), 
values of the internal coordinate. Similar to the pair 
interaction, a permutation of the indices Si in gs^Sji'^) 
results in a rotation in space by tt, so that the functions 
^i^_i(r)and ^_i^i(r) are related by this transformation. 
We obtain these distribution functions at various temper- 
atures and densities in the region where the liquid phase 
is stable by averaging over uncorrelated configurations. 
While ^1,1 (r) is relatively insensitive to temperature, the 
nature of ^i^_i(r) depends strongly on T. We have il- 
lustrated this in Figj7| At low temperatures, when the 
system is dominated by it's potential energy which has 
a 3-fold symmetry, ^i^_i(r) becomes large in magnitude 
compared to ^1,1 (r), long ranged and shows strong direc- 
tionality pointing to the formation of a prominent short 
ranged network in the liquid phase. This network is three 
fold coordinated similar to the HS phase and disappears 
when the temperature is increased. 

In order to further understand our results, we compare 
the pair distribution functions obtained from our simu- 
lations with the results of an approximate integral equa- 
tion theory which we describe below. Unlike a molecular 



fiuid, in our case we need to set up equations for just 
two functions corresponding to like and unlike values of 
Si as in a binary mixture. The correlation functions are 
however direction dependent. In this case, we devise a 
perturbative scheme, where the direction dependent part 
is treated as a "small" perturbation over a set of isotropic 
functions. This allows us to quickly compute distribution 
functions which are in agreement with those obtained 
from simulations at high temperatures, but begins to de- 
viate as the temperature is lowered. Nevertheless, our 
simple scheme is sufficient to show the emergence of the 
short-ranged directional order indicating network forma- 
tion. 

To begin, consider the Ornstein-Zernike equation for a 
binary mixture in Fourier space, viz. 



ha, (30^) =Cc,,/3(k) 



■ X^Ca,j(k)h^,l3{^) 



(8) 



where the indices a, /3 = 1, — 1 denote two species of par- 
ticles, ha, (3 = ga,(3 — 1 are the pair correlation functions, 
Xk are the concentrations of species k (1/2 in our case), 
and Ca,f3 are the direct correlation functions. Note that 
we have suppressed the spatial coordinate for simplicity. 
In order to be useful, EqjS] needs to be supplemented 
(or 'closed') with another equation involving the un- 
known functions ha^/s and Ca^/s- For short-ranged poten- 
tials an approximate closure relation which is known to 
work well is the Percus -Yevick (PY) closure given by. 



Ca,/3(r) =exp[- 



knT 



]{l^yaA^))-yaA^)-l (9) 



with Ha^p = ha,i3 — Ca,(3 the indirect correlation func- 
tion and /7a, /^(r) the pair potential. We shall show that 
within our perturbative scheme, all the correlation func- 
tions may also be decomposed in the same way as Ua,i3 
into isotropic and dependent parts as given in section 
2 A, namely. 
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(10) 



and similarly for y. It is straight forward to show, after 
some algebra, that EqnjS] reduces to. 
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(11) 



where A = \{c'q — Cq) and we have kept only terms up-to 
leading order in Sc. The PY closure, similarly translates 
to, 

Co = exp(^^^j(l + 7/o) - (l + ^o) 




tion data as well as the solutions to the integral equa- 
tions, a strong preference for unlike oriented neighbors 
emerging as the temperature is reduced FigJ8Fc)and(d), 
due to the development of a local network. Since the 
potential energy and hence the 6>-dependent part has a 
significant effect on the system at low temperatures, the 
radial distribution functions obtained from theory and 
simulation however do not match quantitatively in the 
low-temperature limit. More specifically, at tempera- 
tures where our integral equation are valid, the local net- 
work is not strong enough to give rise to thermodynamic 
anomalies which occur for cooler liquids. Better approx- 
imations by keeping orientation dependent terms to all 
orders are needed to yield accurate correlation functions 
in the desired range at the cost of substantially increasing 
the computational complexity. 



FIG. 8: Data obtained from theory (solid curves) and sim- 
ulations(open circles) are plotted together for gi,i{r) for (a) 
T = 3.0, (c) T = 0.5 and gi,-i{r) for (b) T = 3.0, (d) T = 0.5. 



Sc = 



(1 + ^o)[e"^^^o^ cosh{f3SV) - 1] 
-(l + yo)e-^(^o)sinh(/3(^V) 



(12) 



A further, ad- hoc, approximation simplifies the problem 
considerably; the aposteriori justification being given by 
its ability to reproduce some of the more essential fea- 
tures of the correlation functions as compared with the 
output of our simulations. Accordingly, we take SV -> 
in the expression for Cq in Eq|l2| In this case the Eqns. 
[8] and [9] factorize into isotropic and dependent parts. 
We solve the resulting integral equations for the isotropic 
functions self-consistently and then calculate the de- 
pendendent pair correlation functions by iterating EqfT2] 
once. Our results are compared in Figsj7]and[8) In spite 
of the approximations made, we find fair agreement be- 
tween simulations and our integral equation theory espe- 
cially at high temperatures. In FigJ7^(a) and (b) we show 
plots of the radial distribution function gii{r) obtained 
from simulations and solution of the integral equations 
respectively. None of them show any strong directional 
dependence. On the other hand a similar plot of gi-i 
shows a very strong three fold directional dependence 
pointing out the emergence of a local three fold coordi- 
nated, short-ranged network structure. In Fig|8]we plot 
the ^-averaged radial distribution functions ^i,i(r) and 
^1,-1 (^) from theory and simulations both at high and 
low temperature. 

While at high temperature Figj8Fa)and(b) the pair dis- 
tribution function does not show any preference for either 
like or unlike species, one observes in both the simula- 



IV. CONCLUSIONS AND FUTURE 
DIRECTIONS 

We have introduced a model-system in two-dimensions 
to study whether network formation in a system leads 
to polymorphism and thermodynamic anomalies even if 
molecular rotational degrees of freedom are not explicitly 
taken into account. While the connection between direc- 
tional bonding and the existence of multiple states with 
differences in density which causes density anomalies has 
been established quite generally [19], our model shows 
that these properties are robust against rather drastic 
simplifications of the nature of the rotational states. We 
study the equilibrium properties of this system in detail 
and show that, though simple, the model shows many 
features of real liquids like water. 

Our studies should be of direct relevance to the bio- 
logically important case of confined water [26 where ro- 
tations are strongly coupled to translational degrees of 
freedom not unlike the case studied here. Indeed, we 
expect strong network formation in such cases and the 
liquid should show prominent thermodynamic anomalies. 
Calculations in this direction are in progress and will be 
published elsewhere. In future, we aim to use our model 
in more complicated situations to discuss issues such as 
shear fiow and coupling to external fields to further study 
the properties of confined network formers. 
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